The Machine Learning for Social Science Global Subnational Covid Counts Datset

Rex W. Douglass

Introduction

Library Loads

#libraries
library(lubridate)
library(tidyverse)

#devtools::install_github("ropensci/USAboundaries")
#devtools::install_github("ropensci/USAboundariesData")
library(USAboundaries) ; #install.packages('USAboundaries')
data(state_codes)

library(tidyverse)
library(scales)
library(gghighlight)
library(lubridate)
library(R0)  # consider moving all library commands to top -- this one was in a loop below

library(WikidataR)
library(countrycode)

library(usmap) ; # install.packages('usmap')
data(statepop)
#devtools::install_github("ropensci/USAboundaries")
#devtools::install_github("ropensci/USAboundariesData")
library(USAboundaries) ; #install.packages('USAboundaries')
data(state_codes)

library(tidyverse)
library(sf)

library(jsonlite)

#This is too slow it's downloading each
library(GADMTools)
library(strucchange) ; #install.packages('strucchange')
library(tsibble)

library(patchwork)

Data Loading and Cleaning

library(tsibble)

lhs_long <- readRDS("/media/skynet2/905884f0-7546-4273-9061-12a790830beb/rwd_github_private/NESScovid19/data_temp/lhs_long.Rds")
#dim(lhs_long) #187,305

lhs_long_clean <- lhs_long %>%
                  mutate(dataset_gid_geonameid_wikidata_id= paste(dataset,gid,geonameid, wikidata_id, sep="_")  ) %>%
                  mutate(gid_geonameid_wikidata_id= paste(gid,geonameid, wikidata_id, sep="_")  ) %>%
  
                  filter(!is.na(date_asdate)) %>%
                  
                  #First drop anything that's negative. Shouldn't be any negatives.
                  filter(is.na(confirmed) | confirmed>=0) %>%
                  filter(is.na(deaths) | deaths>=0) %>%
                  filter(is.na(tested_people) | tested_people>=0) %>%
                  filter(is.na(tested_samples) | tested_samples>=0) %>%
                  
                  #Then set 0s to NA, we don't trust NAs
                  mutate(confirmed=ifelse(confirmed==0, NA, confirmed)) %>%
                  mutate(deaths=ifelse(deaths==0, NA, deaths)) %>%
                  mutate(tested_people=ifelse(tested_people==0, NA, tested_people)) %>%
                  mutate(tested_samples=ifelse(tested_samples==0, NA, tested_samples)) %>%
                
                  #Then check first differences and drop anything that doesn't weakly increase monotonically
                  group_by(dataset, gid_geonameid_wikidata_id) %>% 
                    arrange(date_asdate) %>%
                    mutate(confirmed_fd=confirmed, deaths_fd=deaths, tested_people_fd=tested_people, tested_samples_fd=tested_samples) %>%
                    tidyr::fill( ends_with("_fd")) %>%
                    mutate_at(vars(ends_with("_fd")), tsibble::difference) %>%
                  ungroup() %>%
                  #About 4k of these observations show a decrease from one time step to the next which suggests either an original error or a joining error
                  #We're going to straight drop those observations as a cleaning step
                  filter( (is.na(confirmed_fd) | confirmed_fd>=0) &  #never allow for a reversal
                          (is.na(deaths_fd) | deaths_fd>=0) &        #never allow for a reversal
                          (is.na(tested_people_fd) | tested_people_fd>=0) &  #never allow for a reversal
                          (is.na(tested_samples_fd) | tested_samples_fd>=0)  #never allow for a reversal
                          ) %>%
  
  
  #we're going to go one step further and require either confirmed or tested to be strictly higher
  #In words, during an episode we don't believe reports mean "no new cases" just no new good reporting
  #metabiota is the worst offender here so restricting it to just that
  filter( !(dataset=="metabiota" & (is.na(confirmed_fd) | confirmed_fd<=0 ))    ) %>% #metabiota is really problematic 90% of what we're doing is trying to account for it
  
  #Also reject any where deaths are greater than confirmed
  filter(is.na(confirmed) | is.na(deaths) | confirmed>=deaths) %>%
  
  #also reject any where confirmed doesn't vary 
  group_by(gid ,  geonameid ,wikidata_id) %>% #chose not to do by dataset because testing datasets might not have confirmed
    filter(var(confirmed, na.rm=T)>0) %>%
  ungroup()
  
#dim(lhs_long_clean) #281,464 #292,580 #299213 #181,563

saveRDS(lhs_long_clean,
        "/media/skynet2/905884f0-7546-4273-9061-12a790830beb/rwd_github_private/NESScovid19/data_temp/lhs_long_cleand.Rds")

Summary Statistics of Our Data

print("Number of observations")

[1] "Number of observations"

dim(lhs_long_clean)

[1] 379074     17

print("Number of locations")

[1] "Number of locations"

lhs_long_clean$wikidata_id %>% unique() %>% length() #3,373

[1] 3696

print("Number of Days")

[1] "Number of Days"

lhs_long_clean$date_asdate %>% unique() %>% length() #113

[1] 125

library(DT) #https://rstudio.github.io/DT/
lhs_long_clean %>% count(dataset) %>% arrange(-n) %>% DT::datatable(options = list(pageLength = 20, autoWidth = TRUE))

#library(gt)
#lhs_long_clean %>% count(dataset) %>% arrange(-n) %>% gt() %>% 
# tab_header(
#    title = md("Location-Days by Dataset")#,
#    #subtitle = "Number of Days with Data by Dataset"
#  ) %>%
#  fmt_missing( columns=everything(), rows=NULL,missing_text = 0)

Spatial Variation in Data Availability

Country Level Data Availability

#gadm36 = st_read("/media/skynet2/905884f0-7546-4273-9061-12a790830beb/rwd_github_private/NESSgadm/data_in/gadm36_gpkg/gadm36.gpkg")
#st_layers("/media/skynet2/905884f0-7546-4273-9061-12a790830beb/rwd_github_private/NESScovid19/data_temp/gadm36_bycountry/gadm36_levels_gpkg/gadm36_levels.gpkg")
gadm36_levels_0 = st_read("/media/skynet2/905884f0-7546-4273-9061-12a790830beb/rwd_github_private/NESScovid19/data_temp/gadm36_bycountry/gadm36_levels_gpkg/gadm36_levels.gpkg", layer="level0")  %>%
                  st_simplify(preserveTopology = FALSE, dTolerance =0.1) #  0.025 this is supposedly broken up by 6 levels and so should have u.s. 
#plot(gadm36_levels_0)
#dim(gadm36_levels_0_sf$sf)
#gadm_plot(gadm36_levels_0_sf)
lhs_long_place_sources <- lhs_long_clean  %>% dplyr::select(gid,   gid_geonameid_wikidata_id) %>% 
                           group_by(gid_geonameid_wikidata_id) %>% count(dataset) %>%
                           group_by(gid_geonameid_wikidata_id) %>%
                           summarise(datasets_n=n()) 

p0 + ggtitle("Covid Count Data Availability at the Country Level")

temp <- gadm36_levels_0 %>%
   left_join(
   lhs_long  %>% count(gid, dataset) %>% rename(GID_0=gid) #there are still dupe geonames wikidata to gid matches that need to be fixed geonameid, wikidata_id, 
) 

Joining, by = "GID_0"
Column `GID_0` joining factor and character vector, coercing into character vector

temp_wide <- temp %>% as.data.frame() %>% dplyr::select(Country=NAME_0, dataset, n) %>% distinct() %>% pivot_wider( id_cols = Country, names_from = dataset,
  values_from = n, values_fill = NULL, values_fn = NULL) 

temp_wide %>% DT::datatable(options = list(pageLength = 20, autoWidth = TRUE))

#install.packages('gt')
#library(gt)
#temp_wide %>% gt() %>% 
# tab_header(
#    title = md("Data Availability by Country"),
#    subtitle = "Number of Days with Data by Dataset"
#  ) %>%
#  fmt_missing( columns=everything(), rows=NULL,missing_text = 0)

State/Province Level Data Availability

gadm36_levels_1 = st_read("/media/skynet2/905884f0-7546-4273-9061-12a790830beb/rwd_github_private/NESScovid19/data_temp/gadm36_bycountry/gadm36_levels_gpkg/gadm36_levels.gpkg", layer="level1")  %>%
  st_simplify(preserveTopology = FALSE, dTolerance =0.1) #  0.025 this is supposedly broken up by 6 levels and so should have u.s. 
p1 <- gadm36_levels_1 %>% 
  left_join(lhs_long_place_sources %>% dplyr::select(gid=gid, datasets_n) %>%  left_join( gadm36_levels_1 %>% as.data.frame() %>% dplyr::select(gid=GID_1, NAME_1) %>% distinct()  )
  ) %>%
  #replace_na(list(datasets_n = 0)) %>% 
  ggplot() + geom_sf(aes(fill = datasets_n)) +
  scale_fill_gradient(low="blue", high="red") +
  theme_bw() 

p1

temp <- gadm36_levels_1 %>%
   left_join(
   lhs_long  %>% count(gid, dataset) %>% rename(GID_1=gid) #there are still dupe geonames wikidata to gid matches that need to be fixed geonameid, wikidata_id, 
) 

Joining, by = "GID_1"
Column `GID_1` joining factor and character vector, coercing into character vector

temp_wide <- temp %>% as.data.frame() %>% dplyr::select(Country=NAME_0, Admin1=NAME_1, dataset, n) %>% distinct() %>% pivot_wider( id_cols = Country:Admin1, names_from = dataset,
  values_from = n, values_fill = NULL, values_fn = NULL) %>% mutate(total = rowSums(.[-c(1,2)], na.rm = T))

temp_wide %>% filter(total>0) %>% DT::datatable(options = list(pageLength = 10, autoWidth = TRUE))

NA

County District Level Data Availability

gadm36_levels_2 = st_read("/media/skynet2/905884f0-7546-4273-9061-12a790830beb/rwd_github_private/NESScovid19/data_temp/gadm36_bycountry/gadm36_levels_gpkg/gadm36_levels.gpkg", layer="level2")  %>%
  st_simplify(preserveTopology = FALSE, dTolerance = 0.001) #  0.025 this is supposedly broken up by 6 levels and so should have u.s. I have to keep shrinking it so misisng goes to zero 
p2 <- gadm36_levels_2 %>% 
  left_join(lhs_long_place_sources %>% dplyr::select(gid=gid, datasets_n) %>%  left_join( gadm36_levels_2 %>% as.data.frame() %>% dplyr::select(gid=GID_2, NAME_2) %>% distinct()  )
  ) %>%
  #replace_na(list(datasets_n = 0)) %>% 
  ggplot() + geom_sf(aes(fill = datasets_n),lwd = 0) +
  scale_fill_gradient(low="blue", high="red") +
  theme_bw() 

p2

temp <- gadm36_levels_2 %>%
   left_join(
   lhs_long  %>% count(gid, dataset) %>% rename(GID_2=gid) #there are still dupe geonames wikidata to gid matches that need to be fixed geonameid, wikidata_id, 
) 

Joining, by = "GID_2"
Column `GID_2` joining factor and character vector, coercing into character vector

temp_wide <- temp %>% as.data.frame() %>% dplyr::select(Country=NAME_0, Admin1=NAME_1, Admin2=NAME_2, dataset, n) %>% distinct() %>% pivot_wider( id_cols = Country:Admin2, names_from = dataset,
  values_from = n, values_fill = NULL, values_fn = NULL) %>% mutate(total = rowSums(.[-c(1,2,3)], na.rm = T))

temp_wide %>% filter(total>0) %>% DT::datatable(options = list(pageLength = 10, autoWidth = TRUE))

It seems your data is too big for client-side DataTables. You may consider server-side processing: https://rstudio.github.io/DT/server.htmlIt seems your data is too big for client-side DataTables. You may consider server-side processing: https://rstudio.github.io/DT/server.html

Data Availability and Interpolation over Time

#test <- lhs_long %>% group_by(dataset, gid, geonameid, wikidata_id, date_asdate) %>% summarize(n=n())

all_na <- function(x) any(!is.na(x)) #https://intellipaat.com/community/12999/remove-columns-from-dataframe-where-all-values-are-na

#bing looks off by a day from the other ones
lhs_wide_qcode <- lhs_long %>% 
                  filter(!is.na(wikidata_id) & !is.na(date_asdate)) %>%
                  group_by(gid, geonameid, wikidata_id, date_asdate) %>%  mutate(confirmed_var=var(confirmed, na.rm=T), deaths_var=var(deaths, na.rm=T)) %>% ungroup() %>%
                  dplyr::select(dataset, gid, geonameid, wikidata_id, date_asdate,confirmed, deaths, tested_samples, tested_people) %>%
                  group_by(dataset, gid, geonameid, wikidata_id, date_asdate) %>%  summarise_if(is.numeric,max,na.rm=T) %>% ungroup() %>% #this is a hack, we should have dupes within datasets
                  pivot_wider(names_from = dataset, values_from = c(confirmed, deaths, tested_samples, tested_people) ) %>% 
                  mutate_if(is.numeric, list(~na_if(abs(.), Inf))) %>% #https://stackoverflow.com/questions/12188509/cleaning-inf-values-from-an-r-dataframe
                  select_if(all_na) 

Here we do the actual interpolation. For each individual location-dataset, we fit a piecewise linear model over time, and then use that to average over day to day noise and interpolate between observations. Because this is in log space, the slope of that line is the day on day percent increase in the count. We finish by fitting a piecewise intercept model to those slopes to get a single daily estimate of how much the count is increasing. Datasets are able to disagree in lots of different ways and we still recover a correct rate of change.

Interpolation of Counts

To each individual dataset-location series, we fit a piecewise linear trend. Each segment is required to have at least 3 observations, but otherwise we make no other constraints. This allows us to flexibly fit linear segments as well as tight curves. We also allow for structural breaks where the series stops, has a structural shift up or down because say a change in reporting rules, but other then otherwise continues as normal. We take the linear trend fit as our best estimate at each point. This averages over day to day noise, e.g. under reporting on the weekend. It also allows us to interpolate between observations where there is missingness in the series but we have reasonable beliefs the true measure continued linearly between the observed points.

dataset_places <- lhs_long_clean %>% dplyr::select(dataset_gid_geonameid_wikidata_id) %>% distinct() %>% pull(dataset_gid_geonameid_wikidata_id)
library(party)
temp_list1 <- list()
for(q in dataset_places  ){
      try({
        temp <- NULL
        temp <- lhs_long_clean %>%
                filter(dataset_gid_geonameid_wikidata_id %in% q) %>%
                arrange(date_asdate) %>%

                mutate(confirmed_log=log(confirmed)) %>%
                mutate(deaths_log=log(deaths)) %>%
                mutate(tested_people_log=log(tested_people)) %>%
                mutate(tested_samples_log=log(tested_samples)) 
          
        if(nrow(temp)==0){next()}
        
        temp <- temp %>% expand(dataset, gid_geonameid_wikidata_id, gid, geonameid,wikidata_id, date_asdate=min(date_asdate):max(date_asdate) %>% as_date() ) %>% 
                 full_join(temp) %>% 
                 mutate(date_asnumeric= as.numeric(date_asdate) ) %>% 
                 arrange(date_asnumeric) %>%
                 mutate(confirmed_nonmissing_count     = (!is.na(confirmed)) %>% cumsum() ) %>%
                 mutate(deaths_nonmissing_count     = (!is.na(deaths)) %>% cumsum() ) %>%
                 mutate(tested_people_nonmissing_count     = (!is.na(tested_people)) %>% cumsum() ) %>%
                 mutate(tested_samples_nonmissing_count     = (!is.na(tested_samples)) %>% cumsum() ) %>%
        
                 arrange(-date_asnumeric) %>%
                 mutate(confirmed_nonmissing_count_reversed     = (!is.na(confirmed)) %>% cumsum() ) %>%
                 mutate(deaths_nonmissing_count_reversed     = (!is.na(deaths)) %>% cumsum() ) %>%
                 mutate(tested_people_nonmissing_count_reversed     = (!is.na(tested_people)) %>% cumsum() ) %>%
                 mutate(tested_samples_nonmissing_count_reversed     = (!is.na(tested_samples)) %>% cumsum() ) %>% 
                 arrange(date_asdate) 

        tryCatch({  
          if(sum(!is.na(temp$confirmed_log))>=3){ #need at least 3 points
            tree_confirmed <- party::mob(confirmed_log ~ 1 + date_asnumeric    | date_asnumeric , control = mob_control(minsplit = 3, alpha=.2), data = temp,  model = linearModel)  # + I(date_asnumeric^2)  
            temp$confirmed_log_y_hat <- predict( tree_confirmed, newdata=temp , type = c("response") )
            temp$confirmed_log_group_number <-  predict( tree_confirmed, newdata=temp , type = c("node"))
            temp <- temp %>% 
                    mutate(confirmed_log_y_hat       =ifelse(confirmed_nonmissing_count>0 & confirmed_nonmissing_count_reversed>0,confirmed_log_y_hat,NA)) %>% #only interpolate within the observed data, never before or after
                    mutate(confirmed_log_group_number=ifelse(confirmed_nonmissing_count>0 & confirmed_nonmissing_count_reversed>0,confirmed_log_group_number,NA)) %>% #only interpolate within the observed data, never before or after

                    group_by(confirmed_log_group_number) %>%
                    arrange(date_asnumeric) %>%
                    mutate(confirmed_log_y_hat_slope=tsibble::difference(confirmed_log_y_hat)) %>% 
                    fill(confirmed_log_y_hat_slope, .direction="up") %>% 
                    ungroup() %>%
                    mutate(confirmed_log_y_hat_percent_change = round((exp(confirmed_log_y_hat_slope)-1)*100,2))  %>%
                    mutate(confirmed_log_y_hat=ifelse(confirmed_log_y_hat<0,NA, confirmed_log_y_hat))
          }
        }, error = function(e) {})
        
        tryCatch({  
          if(sum(!is.na(temp$deaths_log))>=3){ #need at least 3 points
            tree_deaths <- party::mob(deaths_log ~ 1 + date_asnumeric    | date_asnumeric , control = mob_control(minsplit = 3, alpha=.2), data = temp,  model = linearModel)  # + I(date_asnumeric^2)    
            temp$deaths_log_y_hat <- predict( tree_deaths, newdata=temp , type = c("response") )
            temp$deaths_log_group_number <-  predict( tree_deaths, newdata=temp , type = c("node"))
            temp <- temp %>% 
                    mutate(deaths_log_y_hat=ifelse(deaths_nonmissing_count>0 & deaths_nonmissing_count_reversed>0,deaths_log_y_hat,NA)) %>% #only interpolate within the observed data, never before
                    mutate(deaths_log_group_number=ifelse(deaths_nonmissing_count>0 & deaths_nonmissing_count_reversed>0,deaths_log_group_number,NA)) %>% #only interpolate within the observed data, never before

                    group_by(deaths_log_group_number) %>% 
                    arrange(date_asnumeric) %>% 
                    mutate(deaths_log_y_hat_slope=tsibble::difference(deaths_log_y_hat)) %>%
                    fill(deaths_log_y_hat_slope, .direction="up") %>%
                    ungroup() %>%
                    mutate(deaths_log_y_hat_percent_change = round((exp(deaths_log_y_hat_slope)-1)*100,2)) %>%
                    mutate(deaths_log_y_hat=ifelse(deaths_log_y_hat<0,NA, deaths_log_y_hat)) #reject any where the prediction is less than 1 full case
          }
        }, error = function(e) {})
        
        tryCatch({  
          if(sum(!is.na(temp$tested_people_log))>=3){ #need at least 3 points
            tree_tested_people <- party::mob(tested_people_log ~ 1 + date_asnumeric    | date_asnumeric , control = mob_control(minsplit = 3, alpha=.2), data = temp,  model = linearModel)  # + I(date_asnumeric^2)   
            temp$tested_people_log_y_hat <- predict( tree_tested_people, newdata=temp , type = c("response") )
            temp$tested_people_log_group_number <-  predict( tree_tested_people, newdata=temp , type = c("node"))
            temp <- temp %>% 
                       mutate(tested_people_log_y_hat=ifelse(tested_people_nonmissing_count>0 & tested_people_nonmissing_count_reversed>0,tested_people_log_y_hat,NA)) %>% #only interpolate within the observed data, never before
                       mutate(tested_people_log_group_number=ifelse(tested_people_nonmissing_count>0 & tested_people_nonmissing_count_reversed>0,tested_people_log_group_number,NA)) %>% #only interpolate within the observed data, never before
              
                       group_by(tested_people_log_group_number) %>% 
                       arrange(date_asnumeric) %>%
                       mutate(tested_people_log_y_hat_slope=tsibble::difference(tested_people_log_y_hat)) %>% 
                       fill(tested_people_log_y_hat_slope, .direction="up") %>% ungroup() %>%
                       mutate(tested_people_log_y_hat_percent_change = round((exp(tested_people_log_y_hat_slope)-1)*100,2)) %>%
                       mutate(tested_people_log_y_hat=ifelse(tested_people_log_y_hat<0,NA, tested_people_log_y_hat)) #reject any where the prediction is less than 1 full case
          }
        }, error = function(e) {})
        
        temp_list1[[q ]] <- temp
  
      })
}
lhs_long_clean_imputed1 <- bind_rows(temp_list1)
dim(lhs_long_clean_imputed1) #500,786 #bigger because we're interpolating now

#rt1 <- rpart(formula = confirmed_log_y_hat_percent_change ~ date_asnumeric, data=test)
#interpolated = data.frame(date_asnumeric=min(test$date_asnumeric):max(test$date_asnumeric) )
#interpolated$confirmed_log_y_hat_percent_change_y_hat <-  predict(rt1, newdata=interpolated)
#interpolated$date_asdate <- as.Date(interpolated$date_asnumeric)

places <- lhs_long_clean %>% dplyr::select(gid_geonameid_wikidata_id) %>% distinct() %>% pull(gid_geonameid_wikidata_id)

library(rpart)
minsplit=6
minbucket=3
temp_list2 <- list()
for(q in places ){
  print(q)
  temp <- lhs_long_clean_imputed1 %>% 
                 filter(gid_geonameid_wikidata_id %in% q) %>% 
                 arrange(date_asnumeric) 
    
          
  tryCatch({  
    rt1 <- rpart::rpart(formula = confirmed_log_y_hat_percent_change ~ date_asnumeric,
                        data=temp %>% filter(!is.na(confirmed_log)), #to keep from fitting to interpolate slopes that are out of domain and usually wrong, require confirmed to not be missing or na
                        control = rpart.control(minsplit = minsplit, minbucket=minbucket))
    temp$confirmed_log_y_hat_percent_change_y_hat <-  predict(rt1, newdata=temp) #this prediction helpfully tries to impute backwards even before there are any day in any series
    #We need to nuke it if no data exist
    temp <- temp %>%
            group_by(date_asdate) %>%
            mutate( confirmed_nonmissing_count_datasets=sum(confirmed_nonmissing_count>0,na.rm=T) ) %>%
            mutate( confirmed_nonmissing_count_datasets_reversed=sum(confirmed_nonmissing_count_reversed>0,na.rm=T) ) %>%
            mutate(confirmed_log_y_hat_percent_change_y_hat=ifelse(confirmed_nonmissing_count_datasets>0 & confirmed_nonmissing_count_datasets_reversed>0,confirmed_log_y_hat_percent_change_y_hat,NA))  %>%
      ungroup()
    
  }, error = function(e) {})
  
  tryCatch({  
    rt2 <- rpart::rpart(formula = deaths_log_y_hat_percent_change ~ date_asnumeric,
                        data=temp %>% filter(!is.na(deaths_log)),
                        control = rpart.control(minsplit = minsplit, minbucket=minbucket))
    temp$deaths_log_y_hat_percent_change_y_hat <-  predict(rt2, newdata=temp)
    temp <- temp %>%
        group_by(date_asdate) %>%
        mutate( deaths_nonmissing_count_datasets=sum(deaths_nonmissing_count>0,na.rm=T) ) %>%
        mutate( deaths_nonmissing_count_datasets_reversed=sum(deaths_nonmissing_count_reversed>0,na.rm=T) ) %>%
        mutate(deaths_log_y_hat_percent_change_y_hat=ifelse(deaths_nonmissing_count_datasets>0 & deaths_nonmissing_count_datasets_reversed>0,deaths_log_y_hat_percent_change_y_hat,NA)) %>%
      ungroup()
        
  }, error = function(e) {})

  tryCatch({  
    rt3 <- rpart::rpart(formula = tested_people_log_y_hat_percent_change ~ date_asnumeric,
                    data=temp %>% filter(!is.na(tested_people_log)),
                    control = rpart.control(minsplit = minsplit, minbucket=minbucket))
    temp$tested_people_log_y_hat_percent_change_y_hat <-  predict(rt3, newdata=temp, control = rpart.control(minsplit = minsplit, minbucket=minbucket))
    temp <- temp %>%
        group_by(date_asdate) %>%
        mutate( tested_people_nonmissing_count_datasets=sum(tested_people_nonmissing_count>0,na.rm=T) ) %>%
        mutate( tested_people_nonmissing_count_datasets_reversed=sum(tested_people_nonmissing_count_reversed>0,na.rm=T) ) %>%
      
        mutate(tested_people_log_y_hat_percent_change_y_hat=ifelse(tested_people_nonmissing_count_datasets>0 & tested_people_nonmissing_count_datasets_reversed>0,tested_people_log_y_hat_percent_change_y_hat,NA))  %>%
      ungroup()
  }, error = function(e) {})
  
  temp_list2[[q ]] <- temp
}
lhs_long_clean_imputed2 <- bind_rows(temp_list2)

lhs_long_clean_imputed <- lhs_long_clean_imputed2
dim(lhs_long_clean_imputed) #419910 #413475 #408,369 #we lose anywithout q codes here
saveRDS(lhs_long_clean_imputed, "/media/skynet2/905884f0-7546-4273-9061-12a790830beb/rwd_github_private/NESScovid19/data_temp/lhs_long_clean_imputed.Rds")

Examples

China

#China
#"Q148"

test <- lhs_long_clean_imputed %>% 
        mutate(group=paste(dataset,wikidata_id,confirmed_log_y_hat_percent_change)) %>% filter(wikidata_id=="Q148") %>% 
        mutate(group=paste(dataset, wikidata_id,confirmed_log_group_number)) %>%
        arrange(date_asdate) #%>%
        #mutate( confirmed_log_y_hat_percent_change = breakpoints(formula=confirmed_log_y_hat_percent_change ~ 1, data=. ) %>% fitted.values() )

p_interpolation0 <- test %>%
                    ggplot() + 
                       geom_point( aes(x=date_asdate, y=confirmed_log, color=dataset), alpha=.3) +
                       geom_line( aes(x=date_asdate, y=confirmed_log_y_hat, color=dataset, group=group)) + #
                       #geom_line(data=temp_df2, aes(x=date_asdate, y=y_hat_mean), color="black", linetype = "dashed") +
                       theme_bw() +
                       ggtitle("China Q148 Confirmed") + facet_wrap(~dataset   )

Data on China shows variation across datasets for when they start measuring. Only two observe the very beginning of the outbreak. All of the datasets show a discontinuity where reporting standards must have changed. Our method attempts to bridge that gap in 3 of the series but not two others.

p_interpolation0

Italy

#Italy
#"Q38"

test <- lhs_long_clean_imputed %>% mutate(group=paste(dataset,wikidata_id,confirmed_log_y_hat_percent_change)) %>% filter(wikidata_id=="Q38") %>% 
        mutate(group=paste(dataset, wikidata_id,confirmed_log_group_number)) %>%
        arrange(date_asdate) #%>%
        #mutate( confirmed_log_y_hat_percent_change = breakpoints(formula=confirmed_log_y_hat_percent_change ~ 1, data=. ) %>% fitted.values() )

p_interpolation1 <- test %>%
                    ggplot() + 
                       geom_point( aes(x=date_asdate, y=confirmed_log, color=dataset), alpha=.3) +
                       geom_line( aes(x=date_asdate, y=confirmed_log_y_hat, color=dataset, group=group)) + #
                       #geom_line(data=temp_df2, aes(x=date_asdate, y=y_hat_mean), color="black", linetype = "dashed") +
                       theme_bw() +
                       ggtitle("Italy Q38 Confirmed") + facet_wrap(~dataset   )

Italy likewise shows variation in the start of recording across datasets. There’s also variation in the initial sparcity. Wikipedia is strangely missing data in the middle. Our method linearly interpolates into that space in a way in a way we probably don’t want.

p_interpolation1

Italy

#US
#"Q30"

test <- lhs_long_clean_imputed %>% mutate(group=paste(dataset,wikidata_id,confirmed_log_y_hat_percent_change)) %>% filter(wikidata_id=="Q30") %>% 
        mutate(group=paste(dataset, wikidata_id,confirmed_log_group_number)) %>%
        arrange(date_asdate) #%>%
        #mutate( confirmed_log_y_hat_percent_change = breakpoints(formula=confirmed_log_y_hat_percent_change ~ 1, data=. ) %>% fitted.values() )

p_interpolation2 <- test %>%
                    ggplot() + 
                       geom_point( aes(x=date_asdate, y=confirmed_log, color=dataset), alpha=.3) +
                       geom_line( aes(x=date_asdate, y=confirmed_log_y_hat, color=dataset, group=group)) + #
                       #geom_line(data=temp_df2, aes(x=date_asdate, y=y_hat_mean), color="black", linetype = "dashed") +
                       theme_bw() +
                       ggtitle("US Q30 Confirmed") + facet_wrap(~dataset   )

p_interpolation2

New York State

#US
#"Q30"

test <- lhs_long_clean_imputed %>% mutate(group=paste(dataset,wikidata_id,confirmed_log_y_hat_percent_change)) %>% filter(wikidata_id=="Q1384") %>% 
        mutate(group=paste(dataset, wikidata_id,confirmed_log_group_number)) %>%
        arrange(date_asdate) #%>%
        #mutate( confirmed_log_y_hat_percent_change = breakpoints(formula=confirmed_log_y_hat_percent_change ~ 1, data=. ) %>% fitted.values() )

p_interpolation3 <- test %>%
                    ggplot() + 
                       geom_point( aes(x=date_asdate, y=confirmed_log, color=dataset), alpha=.3) +
                       geom_line( aes(x=date_asdate, y=confirmed_log_y_hat, color=dataset, group=group)) + #
                       #geom_line(data=temp_df2, aes(x=date_asdate, y=y_hat_mean), color="black", linetype = "dashed") +
                       theme_bw() +
                       ggtitle("New York State Q1384 Confirmed") + facet_wrap(~dataset   )

p_interpolation3

New York City (by county)

nyt has it just by city, it doesn’t break it down by the boroughs. But because gadm doesn’t do cities we don’t assign it a qcode.

usafacts has it as New York County which is New York County (Manhattan)

#Brooklyn (Q18419)
test <- lhs_long_clean_imputed %>% mutate(group=paste(dataset,wikidata_id,confirmed_log_y_hat_percent_change)) %>% filter(wikidata_id=="Q855974") %>%  #no hits?
        mutate(group=paste(dataset, wikidata_id,confirmed_log_group_number)) %>%
        arrange(date_asdate) #%>%
        #mutate( confirmed_log_y_hat_percent_change = breakpoints(formula=confirmed_log_y_hat_percent_change ~ 1, data=. ) %>% fitted.values() )

p_interpolation4a <- test %>%
                    ggplot() + 
                       geom_point( aes(x=date_asdate, y=confirmed_log, color=dataset), alpha=.3) +
                       geom_line( aes(x=date_asdate, y=confirmed_log_y_hat, color=dataset, group=group)) + #
                       #geom_line(data=temp_df2, aes(x=date_asdate, y=y_hat_mean), color="black", linetype = "dashed") +
                       theme_bw() +
                       ggtitle("Bronx County (Q855974)") + facet_wrap(~dataset   ) + guides( color = FALSE)

#Queens County, New York Q5142559
test <- lhs_long_clean_imputed %>% mutate(group=paste(dataset,wikidata_id,confirmed_log_y_hat_percent_change)) %>% filter(wikidata_id=="Q5142559") %>%  #no hits?
        mutate(group=paste(dataset, wikidata_id,confirmed_log_group_number)) %>%
        arrange(date_asdate) #%>%
        #mutate( confirmed_log_y_hat_percent_change = breakpoints(formula=confirmed_log_y_hat_percent_change ~ 1, data=. ) %>% fitted.values() )

p_interpolation4b <- test %>%
                    ggplot() + 
                       geom_point( aes(x=date_asdate, y=confirmed_log, color=dataset), alpha=.3) +
                       geom_line( aes(x=date_asdate, y=confirmed_log_y_hat, color=dataset, group=group)) + #
                       #geom_line(data=temp_df2, aes(x=date_asdate, y=y_hat_mean), color="black", linetype = "dashed") +
                       theme_bw() +
                       ggtitle("Queens County, New York Q5142559") + facet_wrap(~dataset   ) + guides( color = FALSE)

#Brooklyn is aparently Kings County Q11980692
test <- lhs_long_clean_imputed %>% mutate(group=paste(dataset,wikidata_id,confirmed_log_y_hat_percent_change)) %>% filter(wikidata_id=="Q11980692") %>%  #no hits?
        mutate(group=paste(dataset, wikidata_id,confirmed_log_group_number)) %>%
        arrange(date_asdate) #%>%
        #mutate( confirmed_log_y_hat_percent_change = breakpoints(formula=confirmed_log_y_hat_percent_change ~ 1, data=. ) %>% fitted.values() )

p_interpolation4c <- test %>%
                    ggplot() + 
                       geom_point( aes(x=date_asdate, y=confirmed_log, color=dataset), alpha=.3) +
                       geom_line( aes(x=date_asdate, y=confirmed_log_y_hat, color=dataset, group=group)) + #
                       #geom_line(data=temp_df2, aes(x=date_asdate, y=y_hat_mean), color="black", linetype = "dashed") +
                       theme_bw() +
                       ggtitle("Brooklyn is Kings County Q11980692") + facet_wrap(~dataset   ) + guides( color = FALSE)

# Manhattan is New York County Q500416
test <- lhs_long_clean_imputed %>% mutate(group=paste(dataset,wikidata_id,confirmed_log_y_hat_percent_change)) %>% filter(wikidata_id=="Q500416") %>%  #no hits?
        mutate(group=paste(dataset, wikidata_id,confirmed_log_group_number)) %>%
        arrange(date_asdate) #%>%
        #mutate( confirmed_log_y_hat_percent_change = breakpoints(formula=confirmed_log_y_hat_percent_change ~ 1, data=. ) %>% fitted.values() )

p_interpolation4d <- test %>%
                    ggplot() + 
                       geom_point( aes(x=date_asdate, y=confirmed_log, color=dataset), alpha=.3) +
                       geom_line( aes(x=date_asdate, y=confirmed_log_y_hat, color=dataset, group=group)) + #
                       #geom_line(data=temp_df2, aes(x=date_asdate, y=y_hat_mean), color="black", linetype = "dashed") +
                       theme_bw() +
                       ggtitle("Manhattan is New York County Q500416") + facet_wrap(~dataset   ) + guides( color = FALSE)

#Staten Island is richmond county Q11997784
test <- lhs_long_clean_imputed %>% mutate(group=paste(dataset,wikidata_id,confirmed_log_y_hat_percent_change)) %>% filter(wikidata_id=="Q11997784") %>%  #no hits?
        mutate(group=paste(dataset, wikidata_id,confirmed_log_group_number)) %>%
        arrange(date_asdate) #%>%
        #mutate( confirmed_log_y_hat_percent_change = breakpoints(formula=confirmed_log_y_hat_percent_change ~ 1, data=. ) %>% fitted.values() )

p_interpolation4e <- test %>%
                    ggplot() + 
                       geom_point( aes(x=date_asdate, y=confirmed_log, color=dataset), alpha=.3) +
                       geom_line( aes(x=date_asdate, y=confirmed_log_y_hat, color=dataset, group=group)) + #
                       #geom_line(data=temp_df2, aes(x=date_asdate, y=y_hat_mean), color="black", linetype = "dashed") +
                       theme_bw() +
                       ggtitle("Staten Island  Richmond County Q11997784") + facet_wrap(~dataset   ) + guides( color = FALSE)

Note that CSSE is conflating all of New York City with Manhattan/New York County, Bing is too I think

p_interpolation4a + p_interpolation4b + p_interpolation4c + p_interpolation4d + p_interpolation4e

Interpolation of Rate of Change

Plot Confirmed Curves for Select Locations

library(patchwork)
lhs_long_clean_imputed <- readRDS( "/media/skynet2/905884f0-7546-4273-9061-12a790830beb/rwd_github_private/NESScovid19/data_temp/lhs_long_clean_imputed.Rds")
dim(lhs_long_clean_imputed) #45,609 #337043
#we want the colors to be consistent by dataset too btw

#####
#USA Q30

#bp_temp <- breakpoints(formula=confirmed_log_y_hat_percent_change ~ 1, data=temp_df %>% mutate(group=paste(dataset,wikidata_id,confirmed_log_y_hat_percent_change)) %>% filter(wikidata_id=="Q30"))
#confirmed_log_y_hat_percent_change=fitted.values(bp_temp)

test <- lhs_long_clean_imputed %>%
        mutate(group=paste(dataset,wikidata_id,confirmed_log_group_number)) %>%
        filter(wikidata_id=="Q30") %>% 
        arrange(date_asdate)  %>%
        mutate(date_asnumeric = date_asdate %>% as.numeric() )

#rt1 <- rpart(formula = confirmed_log_y_hat_percent_change ~ date_asnumeric, data=test)
#interpolated = data.frame(date_asnumeric=min(test$date_asnumeric):max(test$date_asnumeric) )
#interpolated$confirmed_log_y_hat_percent_change_y_hat <-  predict(rt1, newdata=interpolated)
#interpolated$date_asdate <- as.Date(interpolated$date_asnumeric)

p0a <- test %>%
        ggplot() + 
           geom_point( aes(x=date_asdate, y=confirmed_log, color=dataset), alpha=.3) +
           geom_line( aes(x=date_asdate, y=confirmed_log_y_hat, color=dataset, group=group)) + #
           #geom_line(data=temp_df2, aes(x=date_asdate, y=y_hat_mean), color="black", linetype = "dashed") +
           theme_bw() +
           ggtitle("USA Q30 Confirmed") # + facet_wrap(~dataset) 

p0b <- test %>%
       ggplot() + 
         geom_point( aes(x=date_asdate , y=confirmed_log_y_hat_percent_change, color=dataset), alpha=.3 ) +
         geom_line(aes(x=date_asdate , y=confirmed_log_y_hat_percent_change_y_hat ), alpha=1, color="black"  ) +
         #geom_smooth( aes(x=date_asdate , y=confirmed_log_y_hat_percent_change), alpha=.3, span=.05) + #, span=.05
         #geom_line( aes(x=date_asdate, y=confirmed_log_y_hat, color=dataset, group=group)) + #
         #geom_line(data=temp_df2, aes(x=date_asdate, y=y_hat_mean), color="black", linetype = "dashed") +
         theme_bw() +
         ggtitle("USA Q30 Confirmed") # + facet_wrap(~dataset)+ ylab("Daily Percent Change Confirmed")

#p0a / p0b



#####
#China Q148

test <- lhs_long_clean_imputed %>% mutate(group=paste(dataset,wikidata_id,confirmed_log_group_number)) %>% filter(wikidata_id=="Q148") %>% 
        mutate(group=paste(dataset, wikidata_id,confirmed_log_group_number)) %>%
        arrange(date_asdate) #%>%

p1a <- test %>% 
ggplot() + 
   geom_point( aes(x=date_asdate, y=confirmed_log, color=dataset), alpha=.3) +
   geom_line( aes(x=date_asdate, y=confirmed_log_y_hat, color=dataset, group=group)) + #
   #geom_line(data=temp_df2, aes(x=date_asdate, y=y_hat_mean), color="black", linetype = "dashed") +
   theme_bw() +
   ggtitle("China Q148 Confirmed")

p1b <- test %>%
       ggplot() + 
         geom_point( aes(x=date_asdate , y=confirmed_log_y_hat_percent_change, color=dataset), alpha=.3) +
         geom_line(aes(x=date_asdate , y=confirmed_log_y_hat_percent_change_y_hat ), alpha=1, color="black") +
         #geom_smooth( aes(x=date_asdate , y=confirmed_log_y_hat_percent_change), alpha=.3, span=.05) + #, span=.05
         #geom_line( aes(x=date_asdate, y=confirmed_log_y_hat, color=dataset, group=group)) + #
         #geom_line(data=temp_df2, aes(x=date_asdate, y=y_hat_mean), color="black", linetype = "dashed") +
         theme_bw() +
   ggtitle("China Q148 Confirmed")+ ylab("Daily Percent Change Confirmed")
   

#p1a / p1b

#India
#"Q668"

test <- lhs_long_clean_imputed %>% mutate(group=paste(dataset,wikidata_id,confirmed_log_y_hat_percent_change)) %>% filter(wikidata_id=="Q668") %>% 
        mutate(group=paste(dataset, wikidata_id,confirmed_log_group_number)) %>%
        arrange(date_asdate) #%>%
        #mutate( confirmed_log_y_hat_percent_change = breakpoints(formula=confirmed_log_y_hat_percent_change ~ 1, data=. ) %>% fitted.values() )


p2a <- test %>%
ggplot() + 
   geom_point( aes(x=date_asdate, y=confirmed_log, color=dataset), alpha=.3) +
   geom_line( aes(x=date_asdate, y=confirmed_log_y_hat, color=dataset, group=group)) + #
   #geom_line(data=temp_df2, aes(x=date_asdate, y=y_hat_mean), color="black", linetype = "dashed") +
   theme_bw() +
   ggtitle("India Q668 Confirmed")

p2b <- test %>%
       ggplot() + 
         geom_point( aes(x=date_asdate , y=confirmed_log_y_hat_percent_change, color=dataset), alpha=.3) +
         geom_line(aes(x=date_asdate , y=confirmed_log_y_hat_percent_change_y_hat ), alpha=1, color="black") +
         #geom_smooth( aes(x=date_asdate , y=confirmed_log_y_hat_percent_change), alpha=.3, span=.05) + #, span=.05
         #geom_line( aes(x=date_asdate, y=confirmed_log_y_hat, color=dataset, group=group)) + #
         #geom_line(data=temp_df2, aes(x=date_asdate, y=y_hat_mean), color="black", linetype = "dashed") +
         theme_bw() +
   ggtitle("India Q668 Confirmed")+ ylab("Daily Percent Change Confirmed")
 
  
#p2a / p2b



#Italy
#"Q38"

test <- lhs_long_clean_imputed %>% mutate(group=paste(dataset,wikidata_id,confirmed_log_y_hat_percent_change)) %>% filter(wikidata_id=="Q38") %>% 
        mutate(group=paste(dataset, wikidata_id,confirmed_log_group_number)) %>%
        arrange(date_asdate) #%>%
        #mutate( confirmed_log_y_hat_percent_change = breakpoints(formula=confirmed_log_y_hat_percent_change ~ 1, data=. ) %>% fitted.values() )


p4a <- test %>%
ggplot() + 
   geom_point( aes(x=date_asdate, y=confirmed_log, color=dataset), alpha=.3) +
   geom_line( aes(x=date_asdate, y=confirmed_log_y_hat, color=dataset, group=group)) + #
   #geom_line(data=temp_df2, aes(x=date_asdate, y=y_hat_mean), color="black", linetype = "dashed") +
   theme_bw() +
   ggtitle("Italy Q38 Confirmed")

p4b <- test %>%
       ggplot() + 
         geom_point( aes(x=date_asdate , y=confirmed_log_y_hat_percent_change, color=dataset), alpha=.3) +
         geom_line(aes(x=date_asdate , y=confirmed_log_y_hat_percent_change_y_hat ), alpha=1, color="black") +
         #geom_smooth( aes(x=date_asdate , y=confirmed_log_y_hat_percent_change), alpha=.3, span=.05) + #, span=.05
         #geom_line( aes(x=date_asdate, y=confirmed_log_y_hat, color=dataset, group=group)) + #
         #geom_line(data=temp_df2, aes(x=date_asdate, y=y_hat_mean), color="black", linetype = "dashed") +
         theme_bw() +
   ggtitle("Italy Q38 Confirmed")+ ylab("Daily Percent Change Confirmed")
   
    
#p4a / p4b



#South Korea
#"Q884"

test <- lhs_long_clean_imputed %>% mutate(group=paste(dataset,wikidata_id,confirmed_log_y_hat_percent_change)) %>% filter(wikidata_id=="Q884") %>% 
        mutate(group=paste(dataset, wikidata_id,confirmed_log_group_number)) %>%
        arrange(date_asdate) #%>%
        #mutate( confirmed_log_y_hat_percent_change = breakpoints(formula=confirmed_log_y_hat_percent_change ~ 1, data=. ) %>% fitted.values() )


p5a <- test %>%
ggplot() + 
   geom_point( aes(x=date_asdate, y=confirmed_log, color=dataset), alpha=.3) +
   geom_line( aes(x=date_asdate, y=confirmed_log_y_hat, color=dataset, group=group)) + #
   #geom_line(data=temp_df2, aes(x=date_asdate, y=y_hat_mean), color="black", linetype = "dashed") +
   theme_bw() +
   ggtitle("South Korea Q884  Confirmed")

p5b <- test %>%
       ggplot() + 
         geom_point( aes(x=date_asdate , y=confirmed_log_y_hat_percent_change, color=dataset), alpha=.3) +
         geom_line(aes(x=date_asdate , y=confirmed_log_y_hat_percent_change_y_hat ), alpha=1, color="black") +
         #geom_smooth( aes(x=date_asdate , y=confirmed_log_y_hat_percent_change), alpha=.3, span=.05) + #, span=.05
         #geom_line( aes(x=date_asdate, y=confirmed_log_y_hat, color=dataset, group=group)) + #
         #geom_line(data=temp_df2, aes(x=date_asdate, y=y_hat_mean), color="black", linetype = "dashed") +
         theme_bw() +
   ggtitle("South Korea Q884 Confirmed") + ylab("Daily Percent Change Confirmed")
   
    
#p5a / p5b

(p0a + p1a + p2a ) / (p0b + p1b + p2b )

( p4a + p5a) / ( p4b + p5b)

test <- lhs_long_clean_imputed %>% mutate(group=paste(dataset,wikidata_id,confirmed_log_y_hat_percent_change)) %>% filter(wikidata_id=="Q1439") %>% 
        mutate(group=paste(dataset, wikidata_id,confirmed_log_group_number)) %>%
        arrange(date_asdate) #%>%
        #mutate( confirmed_log_y_hat_percent_change = breakpoints(formula=confirmed_log_y_hat_percent_change ~ 1, data=. ) %>% fitted.values() )


p0a <- test %>%
      ggplot() + 
         geom_point( aes(x=date_asdate, y=confirmed_log, color=dataset), alpha=.3) +
         geom_line( aes(x=date_asdate, y=confirmed_log_y_hat, color=dataset, group=group)) + #
         #geom_line(data=temp_df2, aes(x=date_asdate, y=y_hat_mean), color="black", linetype = "dashed") +
         theme_bw() +
         ggtitle("Texas Q1439 Confirmed")

p0b <- test %>%
       ggplot() + 
         geom_point( aes(x=date_asdate , y=confirmed_log_y_hat_percent_change, color=dataset), alpha=.3) +
         geom_line(aes(x=date_asdate , y=confirmed_log_y_hat_percent_change_y_hat ), alpha=1, color="black") +
         #geom_smooth( aes(x=date_asdate , y=confirmed_log_y_hat_percent_change), alpha=.3, span=.05) + #, span=.05
         #geom_line( aes(x=date_asdate, y=confirmed_log_y_hat, color=dataset, group=group)) + #
         #geom_line(data=temp_df2, aes(x=date_asdate, y=y_hat_mean), color="black", linetype = "dashed") +
         theme_bw() +
         ggtitle("Texas Q1439 Confirmed") + ylim(0,100)+ ylab("Daily Percent Change Confirmed")

#Bexar
#"Q16861"
test <- lhs_long_clean_imputed %>% mutate(group=paste(dataset,wikidata_id,confirmed_log_y_hat_percent_change)) %>% filter(wikidata_id=="Q16861") %>% 
        mutate(group=paste(dataset, wikidata_id,confirmed_log_group_number)) %>%
        arrange(date_asdate) #%>%
        #mutate( confirmed_log_y_hat_percent_change = breakpoints(formula=confirmed_log_y_hat_percent_change ~ 1, data=. ) %>% fitted.values() )
p1a <- test %>%
      ggplot() + 
         geom_point( aes(x=date_asdate, y=confirmed_log, color=dataset), alpha=.3) +
         geom_line( aes(x=date_asdate, y=confirmed_log_y_hat, color=dataset, group=group)) + #
         #geom_line(data=temp_df2, aes(x=date_asdate, y=y_hat_mean), color="black", linetype = "dashed") +
         theme_bw() +
         ggtitle("Bexar County Q16861 Confirmed")
p1b <- test %>%
       ggplot() + 
         geom_point( aes(x=date_asdate , y=confirmed_log_y_hat_percent_change, color=dataset), alpha=.3) +
         geom_line(aes(x=date_asdate , y=confirmed_log_y_hat_percent_change_y_hat ), alpha=1, color="black") +
         #geom_smooth( aes(x=date_asdate , y=confirmed_log_y_hat_percent_change), alpha=.3, span=.05) + #, span=.05
         #geom_line( aes(x=date_asdate, y=confirmed_log_y_hat, color=dataset, group=group)) + #
         #geom_line(data=temp_df2, aes(x=date_asdate, y=y_hat_mean), color="black", linetype = "dashed") +
         theme_bw() +
   ggtitle("Bexar County Q16861 Confirmed") + ylim(0,100)+ ylab("Daily Percent Change Confirmed")

#California
test <- lhs_long_clean_imputed %>% mutate(group=paste(dataset,wikidata_id,confirmed_log_y_hat_percent_change)) %>% filter(wikidata_id=="Q99") %>% 
        mutate(group=paste(dataset, wikidata_id,confirmed_log_group_number)) %>%
        arrange(date_asdate) #%>%
        #mutate( confirmed_log_y_hat_percent_change = breakpoints(formula=confirmed_log_y_hat_percent_change ~ 1, data=. ) %>% fitted.values() )
p2a <- test %>%
      ggplot() + 
         geom_point( aes(x=date_asdate, y=confirmed_log, color=dataset), alpha=.3) +
         geom_line( aes(x=date_asdate, y=confirmed_log_y_hat, color=dataset, group=group)) + #
         #geom_line(data=temp_df2, aes(x=date_asdate, y=y_hat_mean), color="black", linetype = "dashed") +
         theme_bw() +
         ggtitle("California Q99 Confirmed")
p2b <- test %>%
       ggplot() + 
         geom_point( aes(x=date_asdate , y=confirmed_log_y_hat_percent_change, color=dataset), alpha=.3) +
         geom_line(aes(x=date_asdate , y=confirmed_log_y_hat_percent_change_y_hat ), alpha=1, color="black") +
         #geom_smooth( aes(x=date_asdate , y=confirmed_log_y_hat_percent_change), alpha=.3, span=.05) + #, span=.05
         #geom_line( aes(x=date_asdate, y=confirmed_log_y_hat, color=dataset, group=group)) + #
         #geom_line(data=temp_df2, aes(x=date_asdate, y=y_hat_mean), color="black", linetype = "dashed") +
         theme_bw() +
         ggtitle("California Q99 Confirmed") + ylim(0,100)+ ylab("Daily Percent Change Confirmed")

#San Diego County (Q108143)
test <- lhs_long_clean_imputed %>% mutate(group=paste(dataset,wikidata_id,confirmed_log_y_hat_percent_change)) %>% filter(wikidata_id=="Q108143") %>% 
        mutate(group=paste(dataset, wikidata_id,confirmed_log_group_number)) %>%
        arrange(date_asdate) #%>%
        #mutate( confirmed_log_y_hat_percent_change = breakpoints(formula=confirmed_log_y_hat_percent_change ~ 1, data=. ) %>% fitted.values() )
p3a <- test %>%
      ggplot() + 
         geom_point( aes(x=date_asdate, y=confirmed_log, color=dataset), alpha=.3) +
         geom_line( aes(x=date_asdate, y=confirmed_log_y_hat, color=dataset, group=group)) + #
         #geom_line(data=temp_df2, aes(x=date_asdate, y=y_hat_mean), color="black", linetype = "dashed") +
         theme_bw() +
         ggtitle("San Diego County (Q108143) Confirmed")
p3b <- test %>%
       ggplot() + 
         geom_point( aes(x=date_asdate , y=confirmed_log_y_hat_percent_change, color=dataset), alpha=.3) +
         geom_line(aes(x=date_asdate , y=confirmed_log_y_hat_percent_change_y_hat ), alpha=1, color="black") +
         #geom_smooth( aes(x=date_asdate , y=confirmed_log_y_hat_percent_change), alpha=.3, span=.05) + #, span=.05
         #geom_line( aes(x=date_asdate, y=confirmed_log_y_hat, color=dataset, group=group)) + #
         #geom_line(data=temp_df2, aes(x=date_asdate, y=y_hat_mean), color="black", linetype = "dashed") +
         theme_bw() +
         ggtitle("San Diego County (Q108143) Confirmed") + ylim(0,100)+ ylab("Daily Percent Change Confirmed")

#Florida (Q812)
test <- lhs_long_clean_imputed %>% mutate(group=paste(dataset,wikidata_id,confirmed_log_y_hat_percent_change)) %>% filter(wikidata_id=="Q812") %>% 
        mutate(group=paste(dataset, wikidata_id,confirmed_log_group_number)) %>%
        arrange(date_asdate) #%>%
        #mutate( confirmed_log_y_hat_percent_change = breakpoints(formula=confirmed_log_y_hat_percent_change ~ 1, data=. ) %>% fitted.values() )
p4a <- test %>%
      ggplot() + 
         geom_point( aes(x=date_asdate, y=confirmed_log, color=dataset), alpha=.3) +
         geom_line( aes(x=date_asdate, y=confirmed_log_y_hat, color=dataset, group=group)) + #
         #geom_line(data=temp_df2, aes(x=date_asdate, y=y_hat_mean), color="black", linetype = "dashed") +
         theme_bw() +
         ggtitle("Florida (Q812) Confirmed")
p4b <- test %>%
       ggplot() + 
         geom_point( aes(x=date_asdate , y=confirmed_log_y_hat_percent_change, color=dataset), alpha=.3) +
         geom_line(aes(x=date_asdate , y=confirmed_log_y_hat_percent_change_y_hat ), alpha=1, color="black") +
         #geom_smooth( aes(x=date_asdate , y=confirmed_log_y_hat_percent_change), alpha=.3, span=.05) + #, span=.05
         #geom_line( aes(x=date_asdate, y=confirmed_log_y_hat, color=dataset, group=group)) + #
         #geom_line(data=temp_df2, aes(x=date_asdate, y=y_hat_mean), color="black", linetype = "dashed") +
         theme_bw() +
         ggtitle("Florida (Q812) Confirmed") + ylim(0,100)+ ylab("Daily Percent Change Confirmed")


#St. Lucie County (Q494564)
test <- lhs_long_clean_imputed %>% mutate(group=paste(dataset,wikidata_id,confirmed_log_y_hat_percent_change)) %>% filter(wikidata_id=="Q494564") %>% 
        mutate(group=paste(dataset, wikidata_id,confirmed_log_group_number)) %>%
        arrange(date_asdate) #%>%
        #mutate( confirmed_log_y_hat_percent_change = breakpoints(formula=confirmed_log_y_hat_percent_change ~ 1, data=. ) %>% fitted.values() )
p5a <- test %>%
      ggplot() + 
         geom_point( aes(x=date_asdate, y=confirmed_log, color=dataset), alpha=.3) +
         geom_line( aes(x=date_asdate, y=confirmed_log_y_hat, color=dataset, group=group)) + #
         #geom_line(data=temp_df2, aes(x=date_asdate, y=y_hat_mean), color="black", linetype = "dashed") +
         theme_bw() +
         ggtitle("St. Lucie County (Q494564) Confirmed")
p5b <- test %>%
       ggplot() + 
         geom_point( aes(x=date_asdate , y=confirmed_log_y_hat_percent_change, color=dataset), alpha=.3) +
         geom_line(aes(x=date_asdate , y=confirmed_log_y_hat_percent_change_y_hat ), alpha=1, color="black") +
         #geom_smooth( aes(x=date_asdate , y=confirmed_log_y_hat_percent_change), alpha=.3, span=.05) + #, span=.05
         #geom_line( aes(x=date_asdate, y=confirmed_log_y_hat, color=dataset, group=group)) + #
         #geom_line(data=temp_df2, aes(x=date_asdate, y=y_hat_mean), color="black", linetype = "dashed") +
         theme_bw() +
         ggtitle("St. Lucie County (Q494564) Confirmed") + ylim(0,100)+ ylab("Daily Percent Change Confirmed")

#Not counting slopes that were interpolated out of domain fixed california's plot
( p0a + p1a ) / ( p0b + p1b )

( p2a + p3a ) / ( p2b + p3b )

( p4a + p5a) / ( p4b + p5b)

Plot Death Curves for Select Locations

Show how interpolation works on deaths

test <- lhs_long_clean_imputed %>%
        mutate(group=paste(dataset,wikidata_id,deaths_log_group_number)) %>%
        filter(wikidata_id=="Q30") %>% 
        arrange(date_asdate)  %>%
        mutate(date_asnumeric = date_asdate %>% as.numeric() )

p0a <- test %>%
        ggplot() + 
           geom_point( aes(x=date_asdate, y=deaths_log, color=dataset), alpha=.3) +
           geom_line( aes(x=date_asdate, y=deaths_log_y_hat, color=dataset, group=group)) + #
           #geom_line(data=temp_df2, aes(x=date_asdate, y=y_hat_mean), color="black", linetype = "dashed") +
           theme_bw() +
           ggtitle("USA Q30 Deaths") # + facet_wrap(~dataset) 

p0b <- test %>%
       ggplot() + 
         geom_point( aes(x=date_asdate , y=deaths_log_y_hat_percent_change, color=dataset), alpha=.3 ) +
         geom_line(aes(x=date_asdate , y=deaths_log_y_hat_percent_change_y_hat ), alpha=1, color="black"  ) +
         #geom_smooth( aes(x=date_asdate , y=deaths_log_y_hat_percent_change), alpha=.3, span=.05) + #, span=.05
         #geom_line( aes(x=date_asdate, y=deaths_log_y_hat, color=dataset, group=group)) + #
         #geom_line(data=temp_df2, aes(x=date_asdate, y=y_hat_mean), color="black", linetype = "dashed") +
         theme_bw() +
         ggtitle("USA Q30 Deaths") + ylab("Daily Percent Change Deaths") # + facet_wrap(~dataset)

#p0a / p0b


#####
#China Q148

test <- lhs_long_clean_imputed %>% mutate(group=paste(dataset,wikidata_id,deaths_log_group_number)) %>% filter(wikidata_id=="Q148") %>% 
        mutate(group=paste(dataset, wikidata_id,deaths_log_group_number)) %>%
        arrange(date_asdate) #%>%

p1a <- test %>% 
ggplot() + 
   geom_point( aes(x=date_asdate, y=deaths_log, color=dataset), alpha=.3) +
   geom_line( aes(x=date_asdate, y=deaths_log_y_hat, color=dataset, group=group)) + #
   #geom_line(data=temp_df2, aes(x=date_asdate, y=y_hat_mean), color="black", linetype = "dashed") +
   theme_bw() +
   ggtitle("China Q148 Deaths")

p1b <- test %>%
       ggplot() + 
         geom_point( aes(x=date_asdate , y=deaths_log_y_hat_percent_change, color=dataset), alpha=.3) +
         geom_line(aes(x=date_asdate , y=deaths_log_y_hat_percent_change_y_hat ), alpha=1, color="black") +
         #geom_smooth( aes(x=date_asdate , y=deaths_log_y_hat_percent_change), alpha=.3, span=.05) + #, span=.05
         #geom_line( aes(x=date_asdate, y=deaths_log_y_hat, color=dataset, group=group)) + #
         #geom_line(data=temp_df2, aes(x=date_asdate, y=y_hat_mean), color="black", linetype = "dashed") +
         theme_bw() +
   ggtitle("China Q148 Deaths")+ ylab("Daily Percent Change Deaths")
   

#p1a / p1b

#India
#"Q668"

test <- lhs_long_clean_imputed %>% mutate(group=paste(dataset,wikidata_id,deaths_log_y_hat_percent_change)) %>% filter(wikidata_id=="Q668") %>% 
        mutate(group=paste(dataset, wikidata_id,deaths_log_group_number)) %>%
        arrange(date_asdate) #%>%
        #mutate( deaths_log_y_hat_percent_change = breakpoints(formula=deaths_log_y_hat_percent_change ~ 1, data=. ) %>% fitted.values() )


p2a <- test %>%
ggplot() + 
   geom_point( aes(x=date_asdate, y=deaths_log, color=dataset), alpha=.3) +
   geom_line( aes(x=date_asdate, y=deaths_log_y_hat, color=dataset, group=group)) + #
   #geom_line(data=temp_df2, aes(x=date_asdate, y=y_hat_mean), color="black", linetype = "dashed") +
   theme_bw() +
   ggtitle("India Q668 Deaths")

p2b <- test %>%
       ggplot() + 
         geom_point( aes(x=date_asdate , y=deaths_log_y_hat_percent_change, color=dataset), alpha=.3) +
         geom_line(aes(x=date_asdate , y=deaths_log_y_hat_percent_change_y_hat ), alpha=1, color="black") +
         #geom_smooth( aes(x=date_asdate , y=deaths_log_y_hat_percent_change), alpha=.3, span=.05) + #, span=.05
         #geom_line( aes(x=date_asdate, y=deaths_log_y_hat, color=dataset, group=group)) + #
         #geom_line(data=temp_df2, aes(x=date_asdate, y=y_hat_mean), color="black", linetype = "dashed") +
         theme_bw() +
   ggtitle("India Q668 Deaths")+ ylab("Daily Percent Change Deaths")
 
  
#p2a / p2b

#Bexar
#"Q16861"

test <- lhs_long_clean_imputed %>% mutate(group=paste(dataset,wikidata_id,deaths_log_y_hat_percent_change)) %>% filter(wikidata_id=="Q16861") %>% 
        mutate(group=paste(dataset, wikidata_id,deaths_log_group_number)) %>%
        arrange(date_asdate) #%>%
        #mutate( deaths_log_y_hat_percent_change = breakpoints(formula=deaths_log_y_hat_percent_change ~ 1, data=. ) %>% fitted.values() )


p3a <- test %>%
      ggplot() + 
         geom_point( aes(x=date_asdate, y=deaths_log, color=dataset), alpha=.3) +
         geom_line( aes(x=date_asdate, y=deaths_log_y_hat, color=dataset, group=group)) + #
         #geom_line(data=temp_df2, aes(x=date_asdate, y=y_hat_mean), color="black", linetype = "dashed") +
         theme_bw() +
         ggtitle("Bexar County Q16861 Deaths")

p3b <- test %>%
       ggplot() + 
         geom_point( aes(x=date_asdate , y=deaths_log_y_hat_percent_change, color=dataset), alpha=.3) +
         geom_line(aes(x=date_asdate , y=deaths_log_y_hat_percent_change_y_hat ), alpha=1, color="black") +
         #geom_smooth( aes(x=date_asdate , y=deaths_log_y_hat_percent_change), alpha=.3, span=.05) + #, span=.05
         #geom_line( aes(x=date_asdate, y=deaths_log_y_hat, color=dataset, group=group)) + #
         #geom_line(data=temp_df2, aes(x=date_asdate, y=y_hat_mean), color="black", linetype = "dashed") +
         theme_bw() +
   ggtitle("Bexar County Q16861 Deaths") + ylim(0,100)+ ylab("Daily Percent Change Deaths")
   
#p3a / p3b



#Italy
#"Q38"

test <- lhs_long_clean_imputed %>% mutate(group=paste(dataset,wikidata_id,deaths_log_y_hat_percent_change)) %>% filter(wikidata_id=="Q38") %>% 
        mutate(group=paste(dataset, wikidata_id,deaths_log_group_number)) %>%
        arrange(date_asdate) #%>%
        #mutate( deaths_log_y_hat_percent_change = breakpoints(formula=deaths_log_y_hat_percent_change ~ 1, data=. ) %>% fitted.values() )


p4a <- test %>%
ggplot() + 
   geom_point( aes(x=date_asdate, y=deaths_log, color=dataset), alpha=.3) +
   geom_line( aes(x=date_asdate, y=deaths_log_y_hat, color=dataset, group=group)) + #
   #geom_line(data=temp_df2, aes(x=date_asdate, y=y_hat_mean), color="black", linetype = "dashed") +
   theme_bw() +
   ggtitle("Italy Q38 Deaths")

p4b <- test %>%
       ggplot() + 
         geom_point( aes(x=date_asdate , y=deaths_log_y_hat_percent_change, color=dataset), alpha=.3) +
         geom_line(aes(x=date_asdate , y=deaths_log_y_hat_percent_change_y_hat ), alpha=1, color="black") +
         #geom_smooth( aes(x=date_asdate , y=deaths_log_y_hat_percent_change), alpha=.3, span=.05) + #, span=.05
         #geom_line( aes(x=date_asdate, y=deaths_log_y_hat, color=dataset, group=group)) + #
         #geom_line(data=temp_df2, aes(x=date_asdate, y=y_hat_mean), color="black", linetype = "dashed") +
         theme_bw() +
   ggtitle("Italy Q38 Deaths")+ ylab("Daily Percent Change Deaths")
   
    
#p4a / p4b



#South Korea
#"Q884"

test <- lhs_long_clean_imputed %>% mutate(group=paste(dataset,wikidata_id,deaths_log_y_hat_percent_change)) %>% filter(wikidata_id=="Q884") %>% 
        mutate(group=paste(dataset, wikidata_id,deaths_log_group_number)) %>%
        arrange(date_asdate) #%>%
        #mutate( deaths_log_y_hat_percent_change = breakpoints(formula=deaths_log_y_hat_percent_change ~ 1, data=. ) %>% fitted.values() )


p5a <- test %>%
ggplot() + 
   geom_point( aes(x=date_asdate, y=deaths_log, color=dataset), alpha=.3) +
   geom_line( aes(x=date_asdate, y=deaths_log_y_hat, color=dataset, group=group)) + #
   #geom_line(data=temp_df2, aes(x=date_asdate, y=y_hat_mean), color="black", linetype = "dashed") +
   theme_bw() +
   ggtitle("South Korea Q884  Deaths")

p5b <- test %>%
       ggplot() + 
         geom_point( aes(x=date_asdate , y=deaths_log_y_hat_percent_change, color=dataset), alpha=.3) +
         geom_line(aes(x=date_asdate , y=deaths_log_y_hat_percent_change_y_hat ), alpha=1, color="black") +
         #geom_smooth( aes(x=date_asdate , y=deaths_log_y_hat_percent_change), alpha=.3, span=.05) + #, span=.05
         #geom_line( aes(x=date_asdate, y=deaths_log_y_hat, color=dataset, group=group)) + #
         #geom_line(data=temp_df2, aes(x=date_asdate, y=y_hat_mean), color="black", linetype = "dashed") +
         theme_bw() +
   ggtitle("South Korea Q884 Deaths") + ylab("Daily Percent Change Deaths")
   
    
#p5a / p5b

I need to fix this, USA has a fraction of a death

Note the discontinuity in Chin when they added those 1k deaths to the total

India deaths is broken, we shouldn’t be interpolating a rate before we think the number is greater than 0

(p0a + p1a + p2a ) / (p0b + p1b + p2b )

(p3a + p4a + p5a) / (p3b + p4b + p5b)

Plot Testing Curves for Select Locations

Show how interpolation works on tested

US is an example where we interpolated past our available data. Shouldn’t do that either!

test <- lhs_long_clean_imputed %>%
        mutate(group=paste(dataset,wikidata_id,tested_people_log_group_number)) %>%
        filter(wikidata_id=="Q30") %>% 
        arrange(date_asdate)  %>%
        mutate(date_asnumeric = date_asdate %>% as.numeric() )

p0a <- test %>%
        ggplot() + 
           geom_point( aes(x=date_asdate, y=tested_people_log, color=dataset), alpha=.3) +
           geom_line( aes(x=date_asdate, y=tested_people_log_y_hat, color=dataset, group=group)) + #
           #geom_line(data=temp_df2, aes(x=date_asdate, y=y_hat_mean), color="black", linetype = "dashed") +
           theme_bw() +
           ggtitle("USA Q30 Tested")  #+ facet_wrap(~dataset) 

p0b <- test %>%
       ggplot() + 
         geom_point( aes(x=date_asdate , y=tested_people_log_y_hat_percent_change, color=dataset), alpha=.3 ) +
         geom_line(aes(x=date_asdate , y=tested_people_log_y_hat_percent_change_y_hat ), alpha=1, color="black"  ) +
         #geom_smooth( aes(x=date_asdate , y=tested_people_log_y_hat_percent_change), alpha=.3, span=.05) + #, span=.05
         #geom_line( aes(x=date_asdate, y=tested_people_log_y_hat, color=dataset, group=group)) + #
         #geom_line(data=temp_df2, aes(x=date_asdate, y=y_hat_mean), color="black", linetype = "dashed") +
         theme_bw() +
         ggtitle("USA Q30 tested_people") + ylab("Daily Percent Change tested_people") # + facet_wrap(~dataset)

#p0a / p0b


#####
#China Q148

test <- lhs_long_clean_imputed %>% mutate(group=paste(dataset,wikidata_id,tested_people_log_group_number)) %>% filter(wikidata_id=="Q148") %>% 
        mutate(group=paste(dataset, wikidata_id,tested_people_log_group_number)) %>%
        arrange(date_asdate) #%>%

p1a <- test %>% 
ggplot() + 
   geom_point( aes(x=date_asdate, y=tested_people_log, color=dataset), alpha=.3) +
   geom_line( aes(x=date_asdate, y=tested_people_log_y_hat, color=dataset, group=group)) + #
   #geom_line(data=temp_df2, aes(x=date_asdate, y=y_hat_mean), color="black", linetype = "dashed") +
   theme_bw() +
   ggtitle("China Q148 tested_people")

p1b <- test %>%
       ggplot() + 
         geom_point( aes(x=date_asdate , y=tested_people_log_y_hat_percent_change, color=dataset), alpha=.3) +
         geom_line(aes(x=date_asdate , y=tested_people_log_y_hat_percent_change_y_hat ), alpha=1, color="black") +
         #geom_smooth( aes(x=date_asdate , y=tested_people_log_y_hat_percent_change), alpha=.3, span=.05) + #, span=.05
         #geom_line( aes(x=date_asdate, y=tested_people_log_y_hat, color=dataset, group=group)) + #
         #geom_line(data=temp_df2, aes(x=date_asdate, y=y_hat_mean), color="black", linetype = "dashed") +
         theme_bw() +
   ggtitle("China Q148 tested_people")+ ylab("Daily Percent Change tested_people")
   

#p1a / p1b

#India
#"Q668"

test <- lhs_long_clean_imputed %>% mutate(group=paste(dataset,wikidata_id,tested_people_log_y_hat_percent_change)) %>% filter(wikidata_id=="Q668") %>% 
        mutate(group=paste(dataset, wikidata_id,tested_people_log_group_number)) %>%
        arrange(date_asdate) #%>%
        #mutate( tested_people_log_y_hat_percent_change = breakpoints(formula=tested_people_log_y_hat_percent_change ~ 1, data=. ) %>% fitted.values() )


p2a <- test %>%
ggplot() + 
   geom_point( aes(x=date_asdate, y=tested_people_log, color=dataset), alpha=.3) +
   geom_line( aes(x=date_asdate, y=tested_people_log_y_hat, color=dataset, group=group)) + #
   #geom_line(data=temp_df2, aes(x=date_asdate, y=y_hat_mean), color="black", linetype = "dashed") +
   theme_bw() +
   ggtitle("India Q668 tested_people")

p2b <- test %>%
       ggplot() + 
         geom_point( aes(x=date_asdate , y=tested_people_log_y_hat_percent_change, color=dataset), alpha=.3) +
         geom_line(aes(x=date_asdate , y=tested_people_log_y_hat_percent_change_y_hat ), alpha=1, color="black") +
         #geom_smooth( aes(x=date_asdate , y=tested_people_log_y_hat_percent_change), alpha=.3, span=.05) + #, span=.05
         #geom_line( aes(x=date_asdate, y=tested_people_log_y_hat, color=dataset, group=group)) + #
         #geom_line(data=temp_df2, aes(x=date_asdate, y=y_hat_mean), color="black", linetype = "dashed") +
         theme_bw() +
   ggtitle("India Q668 tested_people")+ ylab("Daily Percent Change tested_people")
 
  
#p2a / p2b

#Bexar
#"Q16861"

test <- lhs_long_clean_imputed %>% mutate(group=paste(dataset,wikidata_id,tested_people_log_y_hat_percent_change)) %>% filter(wikidata_id=="Q16861") %>% 
        mutate(group=paste(dataset, wikidata_id,tested_people_log_group_number)) %>%
        arrange(date_asdate) #%>%
        #mutate( tested_people_log_y_hat_percent_change = breakpoints(formula=tested_people_log_y_hat_percent_change ~ 1, data=. ) %>% fitted.values() )


p3a <- test %>%
      ggplot() + 
         geom_point( aes(x=date_asdate, y=tested_people_log, color=dataset), alpha=.3) +
         geom_line( aes(x=date_asdate, y=tested_people_log_y_hat, color=dataset, group=group)) + #
         #geom_line(data=temp_df2, aes(x=date_asdate, y=y_hat_mean), color="black", linetype = "dashed") +
         theme_bw() +
         ggtitle("Bexar County Q16861 tested_people")

p3b <- test %>%
       ggplot() + 
         geom_point( aes(x=date_asdate , y=tested_people_log_y_hat_percent_change, color=dataset), alpha=.3) +
         geom_line(aes(x=date_asdate , y=tested_people_log_y_hat_percent_change_y_hat ), alpha=1, color="black") +
         #geom_smooth( aes(x=date_asdate , y=tested_people_log_y_hat_percent_change), alpha=.3, span=.05) + #, span=.05
         #geom_line( aes(x=date_asdate, y=tested_people_log_y_hat, color=dataset, group=group)) + #
         #geom_line(data=temp_df2, aes(x=date_asdate, y=y_hat_mean), color="black", linetype = "dashed") +
         theme_bw() +
   ggtitle("Bexar County Q16861 tested_people") + ylim(0,100)+ ylab("Daily Percent Change tested_people")
   
#p3a / p3b



#Italy
#"Q38"

test <- lhs_long_clean_imputed %>% mutate(group=paste(dataset,wikidata_id,tested_people_log_y_hat_percent_change)) %>% filter(wikidata_id=="Q38") %>% 
        mutate(group=paste(dataset, wikidata_id,tested_people_log_group_number)) %>%
        arrange(date_asdate) #%>%
        #mutate( tested_people_log_y_hat_percent_change = breakpoints(formula=tested_people_log_y_hat_percent_change ~ 1, data=. ) %>% fitted.values() )


p4a <- test %>%
ggplot() + 
   geom_point( aes(x=date_asdate, y=tested_people_log, color=dataset), alpha=.3) +
   geom_line( aes(x=date_asdate, y=tested_people_log_y_hat, color=dataset, group=group)) + #
   #geom_line(data=temp_df2, aes(x=date_asdate, y=y_hat_mean), color="black", linetype = "dashed") +
   theme_bw() +
   ggtitle("Italy Q38 tested_people")

p4b <- test %>%
       ggplot() + 
         geom_point( aes(x=date_asdate , y=tested_people_log_y_hat_percent_change, color=dataset), alpha=.3) +
         geom_line(aes(x=date_asdate , y=tested_people_log_y_hat_percent_change_y_hat ), alpha=1, color="black") +
         #geom_smooth( aes(x=date_asdate , y=tested_people_log_y_hat_percent_change), alpha=.3, span=.05) + #, span=.05
         #geom_line( aes(x=date_asdate, y=tested_people_log_y_hat, color=dataset, group=group)) + #
         #geom_line(data=temp_df2, aes(x=date_asdate, y=y_hat_mean), color="black", linetype = "dashed") +
         theme_bw() +
   ggtitle("Italy Q38 tested_people")+ ylab("Daily Percent Change tested_people")
   
    
#p4a / p4b



#South Korea
#"Q884"

test <- lhs_long_clean_imputed %>% mutate(group=paste(dataset,wikidata_id,tested_people_log_y_hat_percent_change)) %>% filter(wikidata_id=="Q884") %>% 
        mutate(group=paste(dataset, wikidata_id,tested_people_log_group_number)) %>%
        arrange(date_asdate) #%>%
        #mutate( tested_people_log_y_hat_percent_change = breakpoints(formula=tested_people_log_y_hat_percent_change ~ 1, data=. ) %>% fitted.values() )


p5a <- test %>%
ggplot() + 
   geom_point( aes(x=date_asdate, y=tested_people_log, color=dataset), alpha=.3) +
   geom_line( aes(x=date_asdate, y=tested_people_log_y_hat, color=dataset, group=group)) + #
   #geom_line(data=temp_df2, aes(x=date_asdate, y=y_hat_mean), color="black", linetype = "dashed") +
   theme_bw() +
   ggtitle("South Korea Q884  tested_people")

p5b <- test %>%
       ggplot() + 
         geom_point( aes(x=date_asdate , y=tested_people_log_y_hat_percent_change, color=dataset), alpha=.3) +
         geom_line(aes(x=date_asdate , y=tested_people_log_y_hat_percent_change_y_hat ), alpha=1, color="black") +
         #geom_smooth( aes(x=date_asdate , y=tested_people_log_y_hat_percent_change), alpha=.3, span=.05) + #, span=.05
         #geom_line( aes(x=date_asdate, y=tested_people_log_y_hat, color=dataset, group=group)) + #
         #geom_line(data=temp_df2, aes(x=date_asdate, y=y_hat_mean), color="black", linetype = "dashed") +
         theme_bw() +
   ggtitle("South Korea Q884 tested_people") + ylab("Daily Percent Change tested_people")
   
    
#p5a / p5b

I need to fix this, USA has a fraction of a death

Note the discontinuity in Chin when they added those 1k deaths to the total

India deaths is broken, we shouldn’t be interpolating a rate before we think the number is greater than 0

(p0a + p1a + p2a ) / (p0b + p1b + p2b )

(p3a + p4a + p5a) / (p3b + p4b + p5b)

Aligning Curves by Takeoff

library(infotheo)
#Once we removed the no variance places my alignment has more mutual information than either first confirmed or 100 confirmed
mutinformation(X=df_slopes %>% dplyr::select(t,t_alligned,t_alligned_firstconfirmed,confirmed_log_y_hat_percent_change_y_hat) %>% discretize(), method="emp") #we outperformed t_alligned_100confirmed but it wasn't avail for most

#df_slopes %>% filter(abs(t_alligned)<20,abs(t_alligned_100confirmed_100)<20) %>% janitor::tabyl(t_alligned, t_alligned_100confirmed_100)
#hist(df_slopes$t_alligned - df_slopes$t_alligned_100confirmed_100, breaks=50)

Plot of curves raw and then curves alligned by takeoff start

The idea here is that completely unconstrained, COVID-19 growth should follow a logistic curve, flat, upswing, constant growth, downswing, and flat again. The assumption of no unconstrained growth no longer holds because the world is reacting, but there’s still a pretty uniformly characteristic upswing across geographic units. We look for this signal of the largest first difference in the daily percent change, call that the takeoff date, and then allign all the time series with that date as 0.

p_t_alligned_firstconfirmed  / p_t_alligned 

What this gives us is a nice normalized dataset where the task is to predict the shape of that distribution. We could try to predict time until the takeoff, the intensity of the growth at the takeoff point, the time until it returns to zero, the area under the curve, or every nook and change.

lhs_long_median %>% ggplot(aes(x=date_asdate, y=confirmed, color=wikidata_id)) + geom_line() + theme(legend.position = "none")

lhs_long_median %>% ggplot(aes(x=date_asdate, y=deaths, color=wikidata_id)) + geom_line() + theme(legend.position = "none")

lhs_long_median %>% ggplot(aes(x=date_asdate, y=deaths, color=wikidata_id)) + geom_line() + theme(legend.position = "none") + scale_y_log10()

library(tsibble)
lhs_long_median %>% arrange(wikidata_id,date_asdate) %>%
  group_by(wikidata_id) %>%
  mutate(deaths_fd=tsibble::difference(deaths)) %>% 
  ungroup() %>%
  filter(deaths_fd>0) %>%
  ggplot(aes(x=date_asdate, y=deaths_fd, color=wikidata_id)) + 
  geom_line() +
  #stat_smooth(method="loess", se = F) +
  theme(legend.position = "none") # + scale_y_log10()

lm1 <- lm(deaths ~  confirmed,data=lhs_long_median %>% filter(CFR<1))
lm2 <- lm(deaths ~ tested_people + confirmed,data=lhs_long_median %>% filter(CFR<1))

library(ggRandomForests); #install.packages('ggRandomForests')
rf1 <- rfsrc(deaths~ confirmed,
                      data=lhs_long_median %>% filter(CFR<1) %>% as.data.frame() )
gg_e <- gg_error(rf1)
gg_v <- gg_variable(rf1)
plot(gg_v, panel=TRUE, se=.95, span=1.2, alpha=.4) 

rf1 <- rfsrc(deaths~ confirmed + tested_people,
                      data=lhs_long_median %>% filter(CFR<1) %>% as.data.frame() )
gg_e <- gg_error(rf1)
gg_v <- gg_variable(rf1)
plot(gg_v, panel=TRUE, se=.95, span=1.2, alpha=.4)


rf2 <- rfsrc(CFR ~  tested_people_log,
                      data=lhs_long_median %>% 
                      mutate(CFR=deaths/confirmed) %>%
                      mutate(tested_people_log=log(tested_people+1)) %>%
                      filter(CFR<1) %>% as.data.frame()
             )
gg_v <- gg_variable(rf2)
plot(gg_v, panel=TRUE, se=.95, span=1.2, alpha=.4)

rf2 <- rfsrc(CFR ~  tested_people_log + confirmed_log,
                      data=lhs_long_median %>% 
                      mutate(CFR=deaths/confirmed) %>%
                      mutate(tested_people_log=log(tested_people+1)) %>%
                      mutate(confirmed_log=log(confirmed+1)) %>%
                      filter(CFR<1) %>% as.data.frame()
             )
gg_v <- gg_variable(rf2)
plot(gg_v, panel=TRUE, se=.95, span=1.2, alpha=.4)


#Throw in time
rf3 <- rfsrc(CFR ~  positive_perc + tested_people_log,
                                    data=lhs_long_median %>% 
                                    mutate(CFR=deaths/confirmed) %>%
                                    mutate(positive_perc=confirmed/tested_people) %>%
                                    mutate(tested_people_log=log(tested_people+1)) %>%
                                    mutate(confirmed_log=log(confirmed+1)) %>%
                                    filter(CFR<1) %>% as.data.frame()
             )
gg_v <- gg_variable(rf3)
plot(gg_v, panel=TRUE, se=.95, span=1.2, alpha=.4)

partial_Boston <- plot.variable(rf3,
partial=TRUE, sorted=FALSE,
show.plots = FALSE )

gg_p <- gg_partial(partial_Boston)
plot(gg_p, panel=TRUE, points = F )


copper_cts <-quantile_pts(lhs_long_median$tested_people_log, groups = 6, intervals = TRUE)
partial_coplot_Boston <- gg_partial_coplot(rf2, xvar="confirmed_log",
                                         groups=rm_grp,
                                         show.plots=FALSE)



summary(lhs_long_median$CFR)

lhs_long_median %>% filter(confirmed>10) %>% filter(CFR<1) %>% pull(CFR) %>% summary() #0.0297297  that's a median CFR of about 3%
lhs_long_median %>% filter(confirmed>10) %>% filter(CFR<1) %>% pull(CFR) %>% hist(breaks=50)

lhs_long_median %>% filter(confirmed>10) %>% filter(CFR<1) %>% ggplot(aes(x=CFR)) + geom_density()

#This code is now depricated in favor of the tree base method above

places <- lhs_long_clean$wikidata_id %>% unique() %>% na.omit() ; length(places)
datasets <- lhs_long_clean$dataset %>% unique() %>% na.omit() ; table(dataset)

temp_list <- list()
for(p in places){
    print(p)
    for(d in datasets){
      temp <- NULL
      temp <- lhs_long_clean %>% 
              arrange(date_asdate) %>% 
              filter(dataset %in% d) %>% 
              filter(wikidata_id %in% p) %>% 
              mutate(confirmed_log=log(confirmed+1)) %>%
              mutate(deaths_log=log(deaths+1)) %>%
              mutate(tested_people_log=log(tested_people+1)) %>%
              mutate(t= as.numeric(date_asdate) - min(as.numeric(date_asdate)) +1  ) %>% #start at 1 to make indexing easier
              mutate(i= row_number()  ) # actually need this
        
      i_original <- temp$i
      
      if( nrow(temp)==0 ) { next} #print("error");
      print(p)
      
      #bp <- breakpoints(confirmed_log ~ 1, data=temp)
      bp <- NULL
      lm1 <- NULL
      y_hat <- NA
      cdf <- NULL
      try({
            #if it fails fall back to just a lm
        lm1 <- lm(confirmed_log ~ 1 + t, data=temp)
    
        cdf <- data.frame(
                          t= 1, 
                          confirmed_log_intercept= coef(lm1)[1],
                          confirmed_log_slope= coef(lm1)[2],
                          confirmed_log_slope_break=0
                          )
        confirmed_log_y_hat=fitted.values(lm1)
      })
      
      try({
        bp <- breakpoints(confirmed_log ~ 1 + t, data=temp, h=3)
        
        cdf <- data.frame(
                          i= c(1,bp$breakpoints), 
                          confirmed_log_intercept= coef(bp)[,1],
                          confirmed_log_slope= coef(bp)[,2],
                          confirmed_log_slope_break=1
                          )
        confirmed_log_y_hat=fitted.values(bp)
      })
      
      try({
        temp2 <- temp
        temp2$confirmed_log_y_hat <- NA
        temp2$confirmed_log_y_hat <- confirmed_log_y_hat
        
        temp2 <- temp2 %>% 
                left_join(cdf) 
        q=paste0(p,"_",d)
        temp2 <- temp2 %>% 
                  expand(dataset, gid,geonameid, wikidata_id,t=min(temp$t):max(temp$t)) %>% #expand this to include all the t 
                  left_join(temp2) %>% 
                  mutate(date_asdate=min(date_asdate, na.rm=T)-1+t) %>% #go back interp date again
                  mutate(confirmed_log_slope_break = confirmed_log_slope_break %>% replace_na(0) %>% cumsum() ) %>%
                  fill(confirmed_log_intercept) %>%
                  fill(confirmed_log_slope) %>%
                  mutate(confirmed_log_y_hat= confirmed_log_intercept + confirmed_log_slope*t ) %>% 
                  mutate(confirmed_log_slope_percent_change = round((exp(confirmed_log_slope)-1)*100,2)) 

        temp_list[[as.character(q)]] <- temp2
      })
      #if( is.na( temp_list[[as.character(q)]]$y_hat) ) {print("error"); break}
    }
}
#"13055"
temp_df <- bind_rows(temp_list)
dim(temp_df) #bigger because we're interpolating now

# install.packages("devtools")
#devtools::install_github("tidyverse/multidplyr")

library(multidplyr)
library(dplyr, warn.conflicts = FALSE)
#cluster <- new_cluster(4)

#need to supress messages
rex_function <- function(x){
  temp=data.frame(confirmed_log_slope=x)
  tryCatch({
    #if there's an NA in x we get an error because breakpoints is one short
    prediction <- breakpoints(formula=confirmed_log_slope ~ 1, data=temp, h=3) %>% fitted.values()
    result=rep(NA, length(x) )
    result[!is.na(x)] <-prediction
    return(result)
    
  }, error=function(e){})
  
  return( rep(NA, length(x) ) )
         
}
#didn't take too long now

#library(multidplyr)
#library(dplyr, warn.conflicts = FALSE)
#cluster <- new_cluster(4)

library(tictoc)
tic()
library(strucchange)
test <- temp_df %>% head(20000)
test3 <- test %>% 
                group_by(wikidata_id) %>% 
                #partition(cluster) %>%
                arrange(date_asdate) %>%
                #mutate( confirmed_log_slope_y_hat = breakpoints(formula=confirmed_log_slope ~ 1, data=. ) %>% fitted.values() ) %>% 
                mutate( confirmed_log_slope_y_hat = rex_function(confirmed_log_slope) ) %>%  #fails on the cluster
                ungroup() #%>%
                #collect()
toc() #48 seconds for 10k #122 seconds for 20k

saveRDS(temp_df, "/media/skynet2/905884f0-7546-4273-9061-12a790830beb/rwd_github_private/NESScovid19/data_temp/lhs_long_interpolated.Rds")